setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
stockdata<-read.table(file="stocks-dat.txt")
colnames(stockdata)<-c("JPMorgan","Citibank","WellsFargo","RoyalDutchShell","ExxonMobil")

## Fit factor model using principal component method with factor rotation
## RC are loadings, h2 are communality, U2 are specific variances
## Com are Hoffman’s index of complexity for each term

library("psych")
PCMfit <- principal(stockdata, nfactors=2, rotate="varimax")
PCMfit

## Possible rotations include: "none", "varimax", "quatimax", "promax", "oblimin", "simplimax", or "cluster"

PCMfitOblimin <- principal(stockdata, nfactors=2, rotate="oblimin")
PCMfitOblimin

IPCMfitOblimin <- fa(stockdata, nfactors=2, rotate="oblimin",fm="wls")
IPCMfitOblimin

## Performs the iterated principal factor method 

IPFfit<-fa(stockdata,nfactors=2, rotate="varimax", fm="wls")
IPFfit

## Maximum likelihood estimator for factor model

MLfit<- factanal(stockdata,2)
print(MLfit, digits=3, cutoff=0, sort=TRUE)
loadings<- MLfit$loadings[,1:2] 
plot(loadings,type="n")
text(loadings,labels=names(stockdata),cex=.7)

MLfitwRotation<-factanal(stockdata,2,rotation="varimax")
MLfitwRotation